Testing Plotting Functions
2023-03-20
TO DO: - add CI’s
Prep
Let’s load a couple random-ish gamlss models and the dataframe they’re built on
## GAMLSS-RS iteration 1: Global Deviance = 1218118
## GAMLSS-RS iteration 2: Global Deviance = 1218033
## GAMLSS-RS iteration 3: Global Deviance = 1218022
## GAMLSS-RS iteration 4: Global Deviance = 1218017
## GAMLSS-RS iteration 5: Global Deviance = 1218017
## GAMLSS-RS iteration 6: Global Deviance = 1218017
## GAMLSS-RS iteration 7: Global Deviance = 1218017
## GAMLSS-RS iteration 8: Global Deviance = 1218017
## GAMLSS-RS iteration 9: Global Deviance = 1218017
## GAMLSS-RS iteration 10: Global Deviance = 1218017
## GAMLSS-RS iteration 11: Global Deviance = 1218017
## GAMLSS-RS iteration 12: Global Deviance = 1218017
## GAMLSS-RS iteration 13: Global Deviance = 1218017Call the script holding the plot functions we want to test:
source("plotting_functions.R")Test
Set-up Functions
These functions are called from the plotting functions.
sim.data() - takes the dataframe you
built your GAMLSS model on and creates a 2 new dfs with simulated data
(male vs female participants) across the age-range, assigning fs_version
and study values to whatever is most common in the original df. Expects
input df to have log_age, fs_version, and
study. Also preps x-axis labels and defines which centiles
you’ll be plotting.
Returns object sim, which is a list of objects.
sim <- sim.data(cn_df)
names(sim)## [1] "ageRange" "dataToPredictM" "dataToPredictF"
## [4] "tickMarks_log" "tickLabels_log" "tickMarks_unscaled"
## [7] "tickLabels_unscaled" "desiredCentiles"centile_predict() - predicts centiles
based on df simulated by sim.data() using the
predictAll() and qGG() functions. Calculates
centiles, 50th centile peak values, and age at peaks separately on male
and female dfs and returns each, as well as an averaged effect across
sexes.
Takes GAMLSS model obj and objects returned from
sim.data(). Returns object pred, which is a
list of objects.
pred <- centile_predict(sGMV.re, sim$dataToPredictM, sim$dataToPredictF, sim$ageRange, sim$desiredCentiles)## new prediction
## new prediction
## new prediction
## new predictionnames(pred)## [1] "fanCentiles" "fanCentiles_M" "fanCentiles_F" "peak"
## [5] "peak_age" "peak_M" "peak_age_M" "peak_F"
## [9] "peak_age_F" "M_mu" "M_sigma" "M_nu"
## [13] "F_mu" "F_sigma" "F_nu"centile_predict.corrected() - same as
centile_predict but corrects out the estimated effects of
fs_version (from mu term) and study (from mu & sigma terms)
parameters.
Expects fs_version to be a fixed effect and
study to be a random effect fit using re()
function!!!
pred.cor <- centile_predict.corrected(sGMV.re, sim$dataToPredictM, sim$dataToPredictF, sim$ageRange, sim$desiredCentiles)## new prediction
## new prediction
## new prediction
## new predictionnames(pred.cor)## [1] "fanCentiles" "fanCentiles_M" "fanCentiles_F" "peak"
## [5] "peak_age" "peak_M" "peak_age_M" "peak_F"
## [9] "peak_age_F" "M_mu" "M_sigma" "M_nu"
## [13] "F_mu" "F_sigma" "F_nu"GGalt.variance - copied from Simon’s Nature paper repo
var <- GGalt.variance(pred$M_mu, pred$M_sigma, pred$M_nu)Plotting
makeCentileFan() - basic centile fan
plotting that averages across sex and predicts on Mode(fs_version) and
Mode(study) of original data the gamlss was modeled on. Expects GAMLSS
model, phenotype being modeled, and the name of the original df.
age_transformed parameter set to TRUE or FALSE
makeCentileFan(GMV.int, "GMV", cn_df, TRUE, "sex")## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
makeCentileFan(GMV.int, "GMV", cn_df, FALSE, "sex")## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
makeCentileFan(sGMV.re, "sGMV", cn_df, TRUE, "sex")## new prediction
## new prediction
## new prediction
## new prediction
makeCentileFan(sGMV.re, "sGMV", cn_df, FALSE, "sex")## new prediction
## new prediction
## new prediction
## new prediction
makeCentileFan_sex_overlay() - same as
makeCentileFan but with separate centile lines for males
and females
makeCentileFan_sex_overlay(GMV.int, "GMV", cn_df, FALSE, "sex")## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
makeCentileFan_sex_overlay(GMV.int, "GMV", cn_df, TRUE, "sex")## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
makeCentileFan_sex_overlay(sGMV.re, "sGMV", cn_df, FALSE, "sex")## new prediction
## new prediction
## new prediction
## new prediction
makeCentileFan_sex_overlay(sGMV.re, "sGMV", cn_df, TRUE, "sex")## new prediction
## new prediction
## new prediction
## new prediction
makeCentileFan.corrected() - centile
fan plot that averages across sexes, correcting for fs_version and study
effects in both centiles and data points plotted.
Expects fs_version to be a fixed effect and
study to be a random effect fit using re()
function!!!
age_transformed parameter set to TRUE or FALSE
makeCentileFan.corrected(sGMV.re, "sGMV", cn_df, TRUE, "sex")## new prediction
## new prediction
## new prediction
## new prediction
#makeCentileFan(sGMV.re, "sGMV", cn_df, TRUE, "sex")
makeCentileFan.corrected(sGMV.re, "sGMV", cn_df, FALSE, "sex")## new prediction
## new prediction
## new prediction
## new prediction
makeCentileFan_sex_overlay.corrected()
- same as makeCentileFan.corrected but with separate
centile lines for males and females
makeCentileFan_sex_overlay.corrected(sGMV.re, "sGMV", cn_df, TRUE, "sex")## new prediction
## new prediction
## new prediction
## new prediction
makeCentileFan_sex_overlay.corrected(sGMV.re, "sGMV", cn_df, FALSE, "sex")## new prediction
## new prediction
## new prediction
## new prediction
plot.gamlss.var() - plots variance from
predicted GAMLSS model separately for males and females.
age_transformed parameter set to TRUE or FALSE
plot.gamlss.var(sGMV.re, "sGMV", cn_df, TRUE)## new prediction
## new prediction
## new prediction
## new prediction
plot.gamlss.var(GMV.int, "GMV", cn_df, FALSE)## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
correct.points() - function used within
site/fs_version corrected plotting to account for these parameter’s
estimate’s effects on each data point. Expects re() to be
used when fitting random effect of site. Can also be called on it’s own.
Returns new df with each subject’s phenotype corrected.
plot_df <- correct.points(sGMV.re, "GMV", cn_df)
summary(plot_df)## GMV sGMV WMV Ventricles
## Min. : 74583 Min. :14132 Min. : 49733 Min. : 811
## 1st Qu.:465826 1st Qu.:52392 1st Qu.: 411303 1st Qu.: 11296
## Median :510770 Median :56277 Median : 451791 Median : 17351
## Mean :519429 Mean :56374 Mean : 453583 Mean : 21542
## 3rd Qu.:565140 3rd Qu.:60497 3rd Qu.: 495940 3rd Qu.: 27358
## Max. :898654 Max. :89344 Max. :1064862 Max. :175416
##
## session participant run study
## 1 :58764 Min. : 1 Min. :1.000 Length:62447
## 2 : 1766 1st Qu.:21374 1st Qu.:1.000 Class :character
## 3 : 1013 Median :41051 Median :1.000 Mode :character
## ses-BL : 111 Mean :40829 Mean :1.025
## ses-v0 : 89 3rd Qu.:61084 3rd Qu.:1.000
## ses-v06: 73 Max. :80991 Max. :8.000
## (Other): 631
## fs_version age_days sex dx country
## Custom : 674 Min. : 291 Female:31646 CN:62447 UK :30081
## FS53 : 3571 1st Qu.: 5393 Male :30801 USA :20734
## FS6_T1 :18219 Median :19273 France : 1816
## FS6_T1T2:39335 Mean :15698 China : 1268
## FSInfant: 648 3rd Qu.:24321 India : 1193
## Max. :36504 Germany: 1036
## (Other): 6319
## site age_yrs log_age TBV
## 11025 :17617 Min. : 0.7967 Min. :2.464 Min. : 136632
## 11027 : 7429 1st Qu.:14.7664 1st Qu.:3.732 1st Qu.: 917810
## 11026 : 3985 Median :52.7659 Median :4.285 Median : 993084
## 1 : 2469 Mean :42.9796 Mean :4.075 Mean : 994554
## MRi-Share: 1341 3rd Qu.:66.5878 3rd Qu.:4.386 3rd Qu.:1072189
## PNC : 1217 Max. :99.9425 Max. :4.562 Max. :1706546
## (Other) :28389
## log_TBV log_GMV log_sGMV log_WMV
## Min. :5.136 Min. :4.873 Min. :4.150 Min. :4.697
## 1st Qu.:5.963 1st Qu.:5.668 1st Qu.:4.719 1st Qu.:5.614
## Median :5.997 Median :5.708 Median :4.750 Median :5.655
## Mean :5.994 Mean :5.710 Mean :4.748 Mean :5.651
## 3rd Qu.:6.030 3rd Qu.:5.752 3rd Qu.:4.782 3rd Qu.:5.695
## Max. :6.232 Max. :5.954 Max. :4.951 Max. :6.027
##
## log_Ventricles site_int pheno_adjust fs_effect
## Min. :2.909 Min. :-0.315826 Min. : 60252 Min. :-0.26998
## 1st Qu.:4.053 1st Qu.:-0.009192 1st Qu.:428767 1st Qu.: 0.07411
## Median :4.239 Median : 0.008215 Median :470435 Median : 0.07865
## Mean :4.250 Mean : 0.001359 Mean :481721 Mean : 0.07444
## 3rd Qu.:4.437 3rd Qu.: 0.008215 3rd Qu.:523826 3rd Qu.: 0.07865
## Max. :5.244 Max. : 0.213378 Max. :840912 Max. : 0.10626
## Experimental
Log-Log(TBV) correction
This function takes log(phenotype) ~ log(TBV) models and plots log(pheno) by age. It predicts on simulated data where for each age, log_TBV is the mean log_TBV (or roughly imputed mean) for that sex at that given age. Hence the “experimental” label, since, as you can see, it’s pretty rough.
## GAMLSS-RS iteration 1: Global Deviance = -324689.2
## GAMLSS-RS iteration 2: Global Deviance = -325130.6
## GAMLSS-RS iteration 3: Global Deviance = -325139.9
## GAMLSS-RS iteration 4: Global Deviance = -325140.3
## GAMLSS-RS iteration 5: Global Deviance = -325140.5
## GAMLSS-RS iteration 6: Global Deviance = -325140.6
## GAMLSS-RS iteration 7: Global Deviance = -325140.7
## GAMLSS-RS iteration 8: Global Deviance = -325140.7
## GAMLSS-RS iteration 9: Global Deviance = -325140.7
## GAMLSS-RS iteration 10: Global Deviance = -325140.8
## GAMLSS-RS iteration 11: Global Deviance = -325140.8
## GAMLSS-RS iteration 12: Global Deviance = -325140.8
## GAMLSS-RS iteration 13: Global Deviance = -325140.8
## GAMLSS-RS iteration 14: Global Deviance = -325140.8
## GAMLSS-RS iteration 15: Global Deviance = -325140.8
## GAMLSS-RS iteration 16: Global Deviance = -325140.7
## GAMLSS-RS iteration 17: Global Deviance = -325140.7
## GAMLSS-RS iteration 18: Global Deviance = -325140.7
## GAMLSS-RS iteration 19: Global Deviance = -325140.7
## GAMLSS-RS iteration 20: Global Deviance = -325140.7
## GAMLSS-RS iteration 21: Global Deviance = -325140.7
## GAMLSS-RS iteration 22: Global Deviance = -325140.6
## GAMLSS-RS iteration 23: Global Deviance = -325140.6
## GAMLSS-RS iteration 24: Global Deviance = -325140.6
## GAMLSS-RS iteration 25: Global Deviance = -325140.6
## GAMLSS-RS iteration 26: Global Deviance = -325140.5
## GAMLSS-RS iteration 27: Global Deviance = -325140.5
## GAMLSS-RS iteration 28: Global Deviance = -325140.5
## GAMLSS-RS iteration 29: Global Deviance = -325140.5
## GAMLSS-RS iteration 30: Global Deviance = -325140.5
## GAMLSS-RS iteration 31: Global Deviance = -325140.4
## GAMLSS-RS iteration 32: Global Deviance = -325140.4
## GAMLSS-RS iteration 33: Global Deviance = -325140.4
## GAMLSS-RS iteration 34: Global Deviance = -325140.4
## GAMLSS-RS iteration 35: Global Deviance = -325140.4
## GAMLSS-RS iteration 36: Global Deviance = -325140.3
## GAMLSS-RS iteration 37: Global Deviance = -325140.3
## GAMLSS-RS iteration 38: Global Deviance = -325140.3
## GAMLSS-RS iteration 39: Global Deviance = -325140.3
## GAMLSS-RS iteration 40: Global Deviance = -325140.3
## GAMLSS-RS iteration 41: Global Deviance = -325140.3
## GAMLSS-RS iteration 42: Global Deviance = -325140.2
## GAMLSS-RS iteration 43: Global Deviance = -325140.2
## GAMLSS-RS iteration 44: Global Deviance = -325140.2
## GAMLSS-RS iteration 45: Global Deviance = -325140.2
## GAMLSS-RS iteration 46: Global Deviance = -325140.2
## GAMLSS-RS iteration 47: Global Deviance = -325140.2
## GAMLSS-RS iteration 48: Global Deviance = -325140.2
## GAMLSS-RS iteration 49: Global Deviance = -325140.1
## GAMLSS-RS iteration 50: Global Deviance = -325140.1
## GAMLSS-RS iteration 51: Global Deviance = -325140.1
## GAMLSS-RS iteration 52: Global Deviance = -325140.1
## GAMLSS-RS iteration 53: Global Deviance = -325140.1
## GAMLSS-RS iteration 54: Global Deviance = -325140.1
## GAMLSS-RS iteration 55: Global Deviance = -325140.1
## GAMLSS-RS iteration 56: Global Deviance = -325140
## GAMLSS-RS iteration 57: Global Deviance = -325140
## GAMLSS-RS iteration 58: Global Deviance = -325140
## GAMLSS-RS iteration 59: Global Deviance = -325140
## GAMLSS-RS iteration 60: Global Deviance = -325140
## GAMLSS-RS iteration 61: Global Deviance = -325140
## GAMLSS-RS iteration 62: Global Deviance = -325140
## GAMLSS-RS iteration 63: Global Deviance = -325140
## GAMLSS-RS iteration 64: Global Deviance = -325140
## GAMLSS-RS iteration 65: Global Deviance = -325140
## GAMLSS-RS iteration 66: Global Deviance = -325139.9
## GAMLSS-RS iteration 67: Global Deviance = -325139.9
## GAMLSS-RS iteration 68: Global Deviance = -325139.9
## GAMLSS-RS iteration 69: Global Deviance = -325139.9
## GAMLSS-RS iteration 70: Global Deviance = -325139.9
## GAMLSS-RS iteration 71: Global Deviance = -325139.9
## GAMLSS-RS iteration 72: Global Deviance = -325139.9
## GAMLSS-RS iteration 73: Global Deviance = -325139.9
## GAMLSS-RS iteration 74: Global Deviance = -325139.9
## GAMLSS-RS iteration 75: Global Deviance = -325139.9
## GAMLSS-RS iteration 76: Global Deviance = -325139.9
## GAMLSS-RS iteration 77: Global Deviance = -325139.9
## GAMLSS-RS iteration 78: Global Deviance = -325139.8
## GAMLSS-RS iteration 79: Global Deviance = -325139.8
## GAMLSS-RS iteration 80: Global Deviance = -325139.8
## GAMLSS-RS iteration 81: Global Deviance = -325139.8
## GAMLSS-RS iteration 82: Global Deviance = -325139.8
## GAMLSS-RS iteration 83: Global Deviance = -325139.8
## GAMLSS-RS iteration 84: Global Deviance = -325139.8
## GAMLSS-RS iteration 85: Global Deviance = -325139.8
## GAMLSS-RS iteration 86: Global Deviance = -325139.8
## GAMLSS-RS iteration 87: Global Deviance = -325139.8
## GAMLSS-RS iteration 88: Global Deviance = -325139.8
## GAMLSS-RS iteration 89: Global Deviance = -325139.8
## GAMLSS-RS iteration 90: Global Deviance = -325139.8
## GAMLSS-RS iteration 91: Global Deviance = -325139.8
## GAMLSS-RS iteration 92: Global Deviance = -325139.8
## GAMLSS-RS iteration 93: Global Deviance = -325139.8
## GAMLSS-RS iteration 94: Global Deviance = -325139.7
## GAMLSS-RS iteration 95: Global Deviance = -325139.7
## GAMLSS-RS iteration 96: Global Deviance = -325139.7
## GAMLSS-RS iteration 97: Global Deviance = -325139.7
## GAMLSS-RS iteration 98: Global Deviance = -325139.7
## GAMLSS-RS iteration 99: Global Deviance = -325139.7
## GAMLSS-RS iteration 100: Global Deviance = -325139.7
## GAMLSS-RS iteration 101: Global Deviance = -325139.7
## GAMLSS-RS iteration 102: Global Deviance = -325139.7
## GAMLSS-RS iteration 103: Global Deviance = -325139.7
## GAMLSS-RS iteration 104: Global Deviance = -325139.7
## GAMLSS-RS iteration 105: Global Deviance = -325139.7
## GAMLSS-RS iteration 106: Global Deviance = -325139.7
## GAMLSS-RS iteration 107: Global Deviance = -325139.7
## GAMLSS-RS iteration 108: Global Deviance = -325139.7
## GAMLSS-RS iteration 109: Global Deviance = -325139.7
## GAMLSS-RS iteration 110: Global Deviance = -325139.7
## GAMLSS-RS iteration 111: Global Deviance = -325139.7
## GAMLSS-RS iteration 112: Global Deviance = -325139.7
## GAMLSS-RS iteration 113: Global Deviance = -325139.7
## GAMLSS-RS iteration 114: Global Deviance = -325139.7
## GAMLSS-RS iteration 115: Global Deviance = -325139.7
## GAMLSS-RS iteration 116: Global Deviance = -325139.7
## GAMLSS-RS iteration 117: Global Deviance = -325139.7
## GAMLSS-RS iteration 118: Global Deviance = -325139.7
## GAMLSS-RS iteration 119: Global Deviance = -325139.6
## GAMLSS-RS iteration 120: Global Deviance = -325139.6
## GAMLSS-RS iteration 121: Global Deviance = -325139.6
## GAMLSS-RS iteration 122: Global Deviance = -325139.6
## GAMLSS-RS iteration 123: Global Deviance = -325139.6
## GAMLSS-RS iteration 124: Global Deviance = -325139.6
## GAMLSS-RS iteration 125: Global Deviance = -325139.6
## GAMLSS-RS iteration 126: Global Deviance = -325139.6
## GAMLSS-RS iteration 127: Global Deviance = -325139.6
## GAMLSS-RS iteration 128: Global Deviance = -325139.6
## GAMLSS-RS iteration 129: Global Deviance = -325139.6
## GAMLSS-RS iteration 130: Global Deviance = -325139.6
## GAMLSS-RS iteration 131: Global Deviance = -325139.6
## GAMLSS-RS iteration 132: Global Deviance = -325139.6
## GAMLSS-RS iteration 133: Global Deviance = -325139.6
## GAMLSS-RS iteration 134: Global Deviance = -325139.6
## GAMLSS-RS iteration 135: Global Deviance = -325139.6
## GAMLSS-RS iteration 136: Global Deviance = -325139.6
## GAMLSS-RS iteration 137: Global Deviance = -325139.6
## GAMLSS-RS iteration 138: Global Deviance = -325139.6
## GAMLSS-RS iteration 139: Global Deviance = -325139.6
## GAMLSS-RS iteration 140: Global Deviance = -325139.6
## GAMLSS-RS iteration 141: Global Deviance = -325139.6
## GAMLSS-RS iteration 142: Global Deviance = -325139.6
## GAMLSS-RS iteration 143: Global Deviance = -325139.6
## GAMLSS-RS iteration 144: Global Deviance = -325139.6
## GAMLSS-RS iteration 145: Global Deviance = -325139.6
## GAMLSS-RS iteration 146: Global Deviance = -325139.6
## GAMLSS-RS iteration 147: Global Deviance = -325139.6
## GAMLSS-RS iteration 148: Global Deviance = -325139.6
## GAMLSS-RS iteration 149: Global Deviance = -325139.6
## GAMLSS-RS iteration 150: Global Deviance = -325139.6
## GAMLSS-RS iteration 151: Global Deviance = -325139.6
## GAMLSS-RS iteration 152: Global Deviance = -325139.6
## GAMLSS-RS iteration 153: Global Deviance = -325139.6
## GAMLSS-RS iteration 154: Global Deviance = -325139.6
## GAMLSS-RS iteration 155: Global Deviance = -325139.6
## GAMLSS-RS iteration 156: Global Deviance = -325139.6
## GAMLSS-RS iteration 157: Global Deviance = -325139.6
## GAMLSS-RS iteration 158: Global Deviance = -325139.6
## GAMLSS-RS iteration 159: Global Deviance = -325139.6
## GAMLSS-RS iteration 160: Global Deviance = -325139.6
## GAMLSS-RS iteration 161: Global Deviance = -325139.6
## GAMLSS-RS iteration 162: Global Deviance = -325139.6
## GAMLSS-RS iteration 163: Global Deviance = -325139.6
## GAMLSS-RS iteration 164: Global Deviance = -325139.6
## GAMLSS-RS iteration 165: Global Deviance = -325139.6
## GAMLSS-RS iteration 166: Global Deviance = -325139.6
## GAMLSS-RS iteration 167: Global Deviance = -325139.6
## GAMLSS-RS iteration 168: Global Deviance = -325139.6
## GAMLSS-RS iteration 169: Global Deviance = -325139.6
## GAMLSS-RS iteration 170: Global Deviance = -325139.6
## GAMLSS-RS iteration 171: Global Deviance = -325139.6
## GAMLSS-RS iteration 172: Global Deviance = -325139.6
## GAMLSS-RS iteration 173: Global Deviance = -325139.5
## GAMLSS-RS iteration 174: Global Deviance = -325139.5makeCentileFan_sex_overlay.logTBV(log_wmv_refit, "log_WMV", cn_df, TRUE, "sex")## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)
## new prediction
## New way of prediction in random() (starting from GAMLSS version 5.0-6)